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1*.  SUPPLEMENTARY  NOTES 

Our  group  uses  theory  and  simulation  as  tools  in  order  to  increase  the 
understanding  of  plasma  instabilities,  heating,  transport,  plasma-wall 
interactions,  and  large  potentials  in  plasmas.  We  also  work  on  the  improvement 
of  simulation  both  theoretically  and  practically. 
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large  potentials  in  plasmas,  bounded  plasmas. 
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This  is  a  brief  progress  report,  covering  our  research  in  general  plasma 
theory  and  simulation,  plasma-wall  physics  theory  and  simulation,  and  code 
development.  Reports  written  in  this  period  are  included  with  this  mailing. 

A  publications  list  plus  abstracts  for  two  major  meetings  are  included. 
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Our  research  for  1989  has  been  widely  reported,  as  given  by  the  listing  following,  of  7  Journal 
Articles,  4  ERL  Reports,  10  Talks,  Conference  Proceedings,  and  8  Talks,  Poster  Papers. 

Abstracts  are  attached  for  most  of  the  talks;  full  papers  for  the  Numerical  Simulation  Conference. 

Sent  along  with  this  Report  are  reprints  of  the  Journal  Articles  and  the  ERL  Reports. 

Our  prior  mode  was  to  publish  Quarterly  Progress  Reports;  these  then  became  Semi-Annual 
Reports,  which  ended  in  1988.  In  1989,  we  began  publishing  Annual  Progress  Reports.  While  QPR’s 
were  excellent  exercises  in  reporting,  they  required  an  immense  effort;  in  today’s  research  climate,  such 
effort  is  not  available. 

We  trust  that  our  reporting  is  still  useful. 
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Publications  for  1989 


Journal  Articles 

W.S.  Lawson,  "Particle  Simulation  of  Bounded  ID  Plasma  Systems,"  J.  Comp.  Physics,  80  (2),  Febru¬ 
ary  1989.  pp.  253-276. 

K.  Theilhaber  and  C.K.  Birdsall,  "Relvin-Helmholtz  Vortex  Formation  and  Particle  Transport  in  a 
Cross-Field  Plasma  Sheath,"  Phys.  Rev.  Lett..  62,  pp.  772-775,  February  13, 1989. 

BJ.  Cohen,  A.B.  Langdon,  D.W.  Hewett  (all  at  LLNL)  aid  RJ.  Procassini  (here),  "Performance  and 
Optimization  of  Direct  Implicit  Particle  Simulation",  J.  Comp.  Pkys.,  81,  pp  151-168,  March  1989 

William  S.  Lawson,  "The  Pierce  Diode  with  an  External  Circuit  I.  Oscillations  About  Nonuniform 
Equilibria,"  Phys.  Fluids  B,  1,  July  1989,  pp  1483-1492. 

William  S.  Lawson,  "The  Pierce  Diode  with  an  External  Circuit  II.  Chaotic  Behavior,”  Phys.  Fluids  B, 
1,  July  1989,  pp.  1493-1501. 

K.  Theilhaber  and  C.K.  Birdsall,  "Relvin-Helmholtz  Vortex  Formation  and  Particle  Transport  in  a 
Cross-Field  Sheath  I:  Transient  Behavior,"  Phys.  Fluids, B,  1,  pp.  2244-2259,  November  1989. 

K.  Theilhaber  and  C.K.  Birdsall,  "Keivin-Helmhoitz  Vortex  Formation  and  Particle  Transport  in  a 
Cross-Field  Sheath  H:  Steady  State,"  Phys.  Fluids  B.  1,  pp2260-2272,  November  1989. 


ERL  Reports 

J.  Verboncoeur,  "ESI  Reference  Manual"  (which  is  distributed  with  our  PC  disk,  which  is  not  included 
here  -  but  free  for  the  asking) 

MJ.  Gerver,  SJE.  Parker,  and  K.  Theilhaber,  "Analytic  Solutions  and  Particle  Simulations  of  Cross-Field 
Plasma  Sheaths."  Memo  No.  UCB/ERL  M89/U4,  August  30,  1989. 

IJ.  Morey,  and  CJC  Birdsall,  "Traveling-Wave-Tube  Simulation:  The  IBC  Code,”  Memo  No. 
UCB/ERL  M89/116,  September  26,  1989. 

SJL  Parker,  and  CJC  Birdsall,  "Numerical  Error  in  Electron  Orbits  with  Large  coc.  At Memo  No. 
UCB/ERL  M89/136,  December  20, 1989. 

Talks,  Conference  Proceedings 

At  Sherwood  Theory  Conference,  San  Antonio,  Texas,  April  3-5,  1989: 

W.S.  Lawson  and  C.K.  Birdsall,  "Simulation  of  RF  driven  plasma  edge." 

SJL  Parker,  C.K.  Birdsall,  A.  Friedman,  S.L.  Ray,  "Bounded  multi-scale  particle  simulation:  The 
sheath  problem." 

At  13th  Coherence  on  Numerical  Simulation  of  Plasmas.  Santa  Fe,  NM,  September  17-20,  1989 

J.P.  Verboncoeur,  V..  Vahedi,  "WinGraphics:  An  Optimized  Windowing  Environment  for  Interac¬ 
tive  Real-Time  Simulations." 

IJ.  Morey,  JJP.  Verboncoeur,  and  V.  Vahedi,  "Bounded  Plasma  Device  Simulation  with  PDW1, 
Including:  External  RLC  Circuit,  DC  and  RF  Drive,  and  Collisional  Processes." 
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IJ.  Morey,  and  C.K.  Birdsall,  "The  Traveling-Wave-Tube  Code  IBC." 

M.V.  Alves,  V.  Vahedi,  and  C.K.  Birdsall,  "PDCl:  One-Dimensional  Radial  Code  for  a  Cylindri¬ 
cal  Plasma  Device  with  an  External  RLC  Circuit” 

SJ5.  Parker,  A.  Friedman,  S.L.  Ray,  and  CJC.  Birdsall,  'Multi-Scale  Particle  Simulation  of 
Bounded  Plasmas." 

S.E.  Parker,  "A  Particle-In-Cell  Method  for  Modeling  Small  Angle  Coulomb  Collisions  in  Plas¬ 
mas.” 

R  J.  Procassini,  C  K.  Birdsall,  B  J.  Cohen,  and  Y.  Matsuda,  'Comparison  of  Particle-In-Cell  and 
Fokker-Planck  Methods  as  Applied  to  the  Modeling  of  Auxiliary-Heated  Mirror  Plasmas." 

Talks,  Poster  Papers 

At  APS/Gaseous  Electronics  Conference  Annual  Meeting,  Palo  Alto,  CA,  October  17-20,  1989 

IJ.  Morey,  V.  Vahedi,  J.P.  Verboncoeur,  and  M.A.  Liebennan,  "Particle  Simulation  Code  for 
Modeling  Processing  Plasmas.” 

M.V.  Alves,  V.  Vahedi,  and  C.K.  Birdsall,  "Cylindrical  Simulations  for  RF  Discharges  and 
Plasma  Immersion  Ion  Implantation.” 

At  APS/Division  of  Plasma  Physics  Annual  Meeting,  Anaheim,  November  13-17,  1989. 

RJ.  Procassini,  C.K.  Birdsall,  and  B.I.  Cohen,  "Particle  Simulations  of  a  Low-Recycling  Divertor 
Scrape-Off  Layer." 

CJC  Birdsall,  RJ.  Procassini,  and  B.I.  Cohen,  "Particle  Simulations  of  a  High-Recycling  Divertor 
Scrape-Off  Layer." 

RJ.  Procassini,  BJ.  Cohen,  Y.  Matsuda,  and  CJC.  Birdsall,  "Modeling  of  Auxiliary-Heated  Mir¬ 
ror  Plasmas:  A  Comparison  of  Particle-In-Cell  and  Fokker-Planck  Methods." 

S.E.  Parker,  and  RJ.  Procassini,  "Large  Space  ami  Tune  Scale  Particle  Simulation  of  Bounded 
Plasmas  with  a  ’Logical  Sheath’." 

M.V.  Alves,  V.  Vahedi,  and  CJC.  Birdsall,  "RF  Plasma  Processes  in  Cylindrical  Models  and 
Plasma  Immersion  Ion  Implantation." 

IJ.  Morey,  V.  Vahedi,  and  J.  Verboncoeur,  "Particle  Simulation  Code  for  Modeling  Processing 
Plasmas." 


Sherwood  Theory  Conference,  San  Antonio, 


Texas,  April  3-5, 


1989 


SW-wood  Theory  Meeting 
Sea  Antonio,  Ttxna 

April  M,  1  m 


Simulation  of  RF  Driven  Plasma  Edge- 
17.  S.  Laamaan  and  C.  JT.  Bird oa£L 

Courut  Institute  <rf  Methcmeticsl  Sdsacro 
New  York  Uairwwtr 
New  York,  NY  10012 

end 

Electroaica  Re— rek  Lehomtory 
Uoirernty  of  CaiHoroie 
Berkeley,  CA  MT20 

Electroetmiic  computer  AmUIobi  of  e  plasma  edge  drives  by  aa  RF  electric 

field  ;n  e  f-p  i  * — rh-tTn  (FrorJij- ekaetda)  effl  be  prienfed.  Result* 

•bow  that  both  the  ioo  flux  and  the  tea  imped  euipe  mrm  raaghly  proportion*] 
to  the  Applied  power.  Both  the  flax  aad  average  imped  energy  are  strong  func¬ 
tions  at  the  drive  frequency  Cor  the  wide  range  at  frequencies  «— H,  but  may  vary 
little  aver  the  narrower  range  at  frequences  need  in  czpertssmsU.  ton  heeling  oc¬ 
curs  over  the  entire  rimulaiion  region  implying  that  wuvee  (loo  Bcmstan  warm) 
ere  directly  excited  at  the  aatenne.  Direct  obs  avail  on  ot  throe  waves  will  be 
attempted,  and  the  reruita  will  be  pre— ted. 


*Thi*  research  wee  supported  by  U.  3.  Department  of  Energy  grants  no.  DE- 
FG02-49ER33223  end  no.  DE-FC01-MER3322Q. 


Bounded  Multi-Scale  Partlde  Simulation:  The  Sheath  Problem 
X  A  Mar  mdC.tC 

HectrvdMU  lUamrK A  IshiUiq,  7n<v.  of  Ct/i/ercs,  StrimUy 
A.  frimimmn  mU  S.  L.  SUt j 

Laarmm  £iwm«t  ifaOoiaal  L atsufnn,  Umm.  of  California.  Lima rwet 

W*  in  roiag  the  nald-scele  [l]  taduuqa*  to  ttuxid  bounded.  r>  items.  Certain  boeadsd 
ryiums  in  &  airmail y  raited  (or  the  melci- teal*  method  becuue  the  boaadtry  Uyvr  (or  iheach) 
that  form!  at  the  wait  which  ii  uaaBy  a  short  tpatul(~  Ao.)  and  dae  scale  (~^)  scrartare. 
can  usxuA.AAdy  effect  the  balk  plums  bekivtor.  One  goal  U  to  udenaad  the  intencuee 
becwvce  the  bulk  plasma  and  the  ihtaih.  If  the  relevant  thort  time  scale  physics  is  local  to  a 
few  (or  one)  spadal  regions,  then  one  oa  take  advantage  of  this  by  advaaciag  partidro  with 
variable  At  depending  oe  position,  hence  rvdedng  competing  time.  The  unmsgnedzed  theath 
problem,  is  ««ch  a  cum. 

The  model  U  a  one  dimensional  bounded  dah  with  kinetic  ions  sad  electrons.  W«  ,«->rT  *,(]} 
a  coUoioulesi  and  onmagnedted  lyiteta  for  tuapUdty.  The  dght  boaxdary  is  a  condemn;  wall 
that  lbvorbs  aA  partidea  that  come  La  contact  with  it.  The  left  boundary  U  a  symmetry  point, 
where  che  partidea  are  reflected.  W«  allow  a  rpectfled  initial  distribution:  f[x,  y,t  m  fl),  sad  a 
dlscnbaced  source:  s(e.e.l). 

To  '.est  the  eameria  of  both  the  multi-scale  method  and  boundary  conditions  w  m  usieg 
the  following  test  problem:  a  cutoff  hlarweQiaa  diatnbuuon  for  the  electrons  and  Ixed  (or 
Inihijteiy  aaaatro)  Ions.  The  system  has  aa  aaalt-nc  solution,  so  the  run  may  be  started  from 
equilibrium.  This  give*  us  a  benchmark  to  compare  against  aad  inti  she  fast  time  scale  electron 
iknnrh  dynamics,  lei  da!  resniU  ualag  roriahle  At  (*uh  ratio  of  amaJUst  to  largest  of  1:32)  will 
be  pveu.  la  the  As  cure,  we  tatend  to  nee  the  more  general  model  to  study  time  dependent 
sheath  problems,  such  ns  a  shock  front  interaction  with  the  boundary. 

[ij  "P^rude-la-Cnll  Plasma  Simulation  with  a  Wide  Range  of  Space  and  Time  Scales*. 

A_  Friedman,  5.L.  Ray,  C.K.  Biriiall.  and  S.£.  Parker,  Proc.  12th  CoeC  ou  .Hummed 
Simulation  of  Plasmas.  APS,  San  Fras cisco,  Sept.  1388. 
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APS/Division  of  Plasma  Physics  Annual.  Meeting,  Anaheim, CA,  November  13-17,  1989 


Particle  Simulations  of  a  Low-Recycling  Diver¬ 
tor  Scrape-Off  Layer  7  R.  Phocajsini,  C.  BlltOSAlL,  Uniter- 
•ill  of  Californio.  Berkeley,  B.  ContN,  Laurence  Livermore  National 
L oho  ratory—1 The  effect  of  Coulomb  collisionallty  on  the  transport  of  par¬ 
ticle,  and  energy  to  a  Iow-reryding  divertor  plate  through  the  icrape-aff 
layer  (SOL)  in  a  tokamak  la  studied  using  a  fully-kinetic,  self-consistent 
particle-la-cell  (PIC)  model.  The  two  specie,  (electron  aad  ion)  PIC 
code  feature,  a  Moate  Carlo  binary-particle  collision  model  which  con¬ 
serve,  the  momentum  and  kinetic  energy  of  the  interacting  particle,. 
The  full-range  of  Coulomb  colliiiontl  processes  (e-«,  t-i  and  i-i)  are  in¬ 
cluded.  The  dependence  of  the  preiheath  and  collector  lheath  potential 
drop,,  plasma  tempecature,  flow  velocity  and  par-dial  heat  flux  on  the 
col!i,ionality  I,  discuiied.  Tlie  “colliilonal"  preaheath  drop,  cl-ctron 
temperature  at  the  plate,  and  collector  ,healh  drop  Increase  with  col- 
lisonality  The  electron  temperature  at  the  plate  ii  much  smaller  than 
the  source  temperature,  due  to  the  lo,i  of  electron*  over  the  collector 
•heath  potential  barrier.  The  electroo  heat  flux,  which  i,  the  main  chan¬ 
nel  for  »n-rry  tran«porl  along  the  Held  line,  in  the  SOL,  It  reduced  rela¬ 
tive  to  the  Treu-.treaming’  heat  flux  y,  2  n.m.ItT*,  for  Ante  (non-mo) 
value,  of  the  mean-free  path  X  a  0,,/u^u. 

*Tlu,  work  supported  by  the  1/  3.  OOE  under  contract  W-MOl-Enq-  ID. 


Particle  Simulations  of  a  High- Recycling  Diver¬ 
tor  Scrape-Off  Layer  7  C.  Bixosall,  R.  PxocASStNi,  Uniter- 
ritj  of  California,  Berkeley,  B.  CotUN,  Lawrence  Lieermore  National 
Laioretory—' The  effect  of  neutral/charged  particle  interactions  on  the 
transport  of  particles  and  energy  to  a  high-recycling  divertor  plate 
throurh  the  scrape-off  layer  (SOL)  in  a  tokamak  it  studied  uting  a  (ully- 
kinetic,  self-cool  it  teat  p&rticle-in-ceU  (PIC)  modeL  In  addition  to  a 
Monte  Carlo  binary-particle  Coulomb-collision  model,  the  basic  elec¬ 
tro,  u  c  particle  code  also  simulates  the  charge-exchange  and  impact- 
lonixaiton  processes  which  occur  between  plasma  particle,  and  recycled 
neutral  particle*  in  the  vicinity  of  the  divertor  plate.  (The  neutral  par¬ 
ticle  d'nsity  proflte  it  specified  by  the  user,  hence  the  neutral  particles 
are  nof  hsndled  In  a  self-consisttnt  manner).  Tbtit  atomic  physics  mod¬ 
els  mry  also  be  used  at  a  source  or  sputtered  impurity  particles,  which 
can  thro  be  followed  by  the  code  as  a  second  ionic  species.  The  de¬ 
pend; -ice  of  tb*  presheath  and  collector  sheath  potential  drops,  plasma 
temperature,  flow  velocity  and  parallel  heat  flux  on  the  ratio  or  the 
char?-  exehange/ionitation  rate  to  the  Coulomb  eotlision  frequency  is 
dlscuswid.  Comparison  it  made  to  other  high-recycling  divertor  models. 

•This  »ork  supported  by  the  U.S.  DOE  under  contract  W-7403-Eng-4S. 


APS/Division  of  Plasma  Physics  Annual  Meeting 
Anaheim,  CA,  November  13-17,  1989 


Modeling  of  Auxiliary-Heated  Mirror  Plasmas: 
A  Comparison  of  Partide-Io-Csll  and  Fokksr-Planclc 
Methods  ‘  J.  PftocASsnn,  B.  Co  him,  Y.  Matjuda,  Uwmet 
Lhermm  National  Laboratory-,  C.  BOUMALL,  Unnanity  of  Califor- 
s is,  Berkeley  —The  transport  and  confinement  of  charted  particles  in 
an  aiudliary-heeced  mirror  plasma  is  modeled  via  a  bounce-averaged 
Fokker- Planck  (F-P)  coda  and  a  direct- imp  licit  particle- in-cell  (PIC) 
code.  The  teat  caee  studied  it  that  at  a  taadam  minor  and  plug  plasma 
which  is  haatad  by  the  injection  of  fundamental  and  second  harmonic 
ECRfl  wave  energy.  This  test  cate  detannines  the  confinement  and 
transport  of  electrons  only,  with  the  ions  modeled  at  scattering  sitae. 
Both  electron-election  and  electrau-iou  collisions  an  included.  The  mag¬ 
netic  and  electrostatic  fields  an  prescribed  quantities.  Each  code  em¬ 
ploys  a  relativistic  description  of  the  electron  dynamics  in  ona  spatial 
dimension.  The  modeling  comparison  is  divided  into  three  sections:  i) 
Benchmarking  the  physics  results  from  the  PIC  code  against  those  from 
the  P-P  code;  a)  A  computational  cost  analysis  of  each  code;  Hi)  A 
discussion  of  the  advantages  and  disadvantages  of  each  code,  including 
a  comparison  of  the  effort  involved  in  setting  up  each  input  deck. 

This  work  supported  by  the  C  ".  DOE  under  contract  W-740 5- Eng-48. 


Particle  Simulation  Coda  for  Modeling  Process¬ 
ing  Plasmas  T  U.  Mow,  V.  Vaicdi  AMD  J.  Vcuoncoidk, 
ffnwersity  of  California,  Berkeley— The  bounded  plasma  particle  simu¬ 
lation  coda  PDW11  has  beta  l  to  ran  interactively  on  a  PC.  A 

collision  model  baa  beta  added  which  can  be  used  to  simulate  elastic, 
esmtarion,  ionisation  and  charge  exchange  collisions.  The  coda  already 
had  models  at  external  circuits  and  plasma  sources,  so  a  number  of  dif- 
ftrmt  processing  pUisui  cut  sow  bo  This  includes  RP  and 

DC  discharges  and  plasma  immersion  ion  implantation  devices.  RF  dis¬ 
charge  Simula  cion,  have  shown  rectification  of  the  plasma  potential,  ions 
responding  to  the  average  potential,  shaath  heating  and  joule  heating. 
Both  voltage  and  current  driven  BP  discharges  have  been  examined, 
and  differences  have  been  observed  between  the  harmonic  content  of  the 
potential.  Ion  velocity  distributions  at  the  boundaries  exhibit  features 
due  to  charge  exchange  and  ionization  within  the  sheath  regions.  The 
cade  has  excellent  graphics  aad  the  many  diagnostics  are  displayed  in  a 
window  format,  accessed  axing  a  menu  system. 

This  work  supported  in  part  by  DOE  contract  DE-FG03-36ER53220 

and  ONR  contract  N00014-85-K-0809 

*WF.  Lawson,  J.  Comp.  Phys.  80,  253  (1989). 


Largo  space  anti  time  scale  particle  simulation  of 
bounded  plasmas  with  n  “Logical  Sheath”?  S.  E.  Paikci 
and  R.  PtOCAMtNl,  ffwusrssty  of  California,  Berkeley  "  When  study¬ 
ing  plssma  flow  to  a  conducting,  absorbing  plate,  the  sheath  can  affect 
the  relevant  physics.  The  sheath  is  a  small  space  (~  Ao.)  and  time 
(~  w~’)  scale  structure.  Hence,  to  inriade  boundary  effect!  in  a  particle 
omulatioa  ordinarily  one  mat  resolve  these  small  time  end  space  ecalee. 
Fortunately,  the  Direct  Implicit  method  has  relaxed  the  uyAt  restraint 
for  stability,  but  is  not  accurate  at  large  w^At  for  bounded  simulations 
in  which  awe  wants  to  retain  the  sheath  dynamics.  We  propose  an  alter¬ 
native  to  resolving  the  sheath  by  implementing  what  we  call  the  ‘logical 
shaath*.  We  absorb  ions  sad  reflect  electrons  sequentially  to  maintain 
mo  net  entrant  to  the  sralL  The  sheath  drop  it  retained  through  the 
cutoff  velocity  of  the  electrons.  Small  time  and  space  scale  resolution 
I*  not  necessary.  We  have  already  used  the  logical  sheath  in  explicit 
simulations  of  a  Q-machine.  We  plan  to  show  implicit  results  with  large 
Up, At.  By  implementing  the  logical  sheath,  we  plan  to  study  kinetic 
effects  tat  the  divertor  region  of  a  tokamak. 

This  work  supported  by  DoE  contract  .Vo.  DE-FG03-86ER53220  and 
OffR  contract  N0001445-K-0809. 


RF  Plasma  Processes  in  Cylindrical  Models  and 
Plasma  Immersion  Ion  Implantation  "  M.  V.  Alves  !  V.  Va- 
■IDI  AMS  C.  K.  BtnoSALL,  ERL,  University  of  California,  Berkeley— 
Spherical  aad  cylindrical  maay-partide  models  are  being  used  to 
eimulato  RF  dischargee  in  which  the  RP  powered  and  the  pounded 
electrodes  have  different  areas.  Then  also  model  plasma  immersion 
km  implantation,  where  a  target  is  pulsed  negatively.  FOCI,  is  a  ID 
radial,  electrostatic  code  rimolating  a  plasma  contained  between  concen¬ 
tric  cylinders  coupled  to  an  external  RLC  circuit  aad  an  BP  source  Col¬ 
lisions  with  neutral  particlei  an  included  as  in  PDWl1  so  that  weakly 
ionized  discharges  can  be  simulated.  The  self-bias  voltage  (producing 
the  ion  bombarding  energy)  at  the  powered  electrode  is  measured  aad 
compared  with  theory.1  Pulsed  plasma  immersion  ion  implantation  will 
also  be  described. 

This  work  supported  is  part  by  DOE  contract  DE-FG03-S6ER53220 
aad  OffR  contract  X0001445-EC-0809 
'On  leave  from  DfPE  -  S.  3.  Campos  -  SP  ■  Brazil 
*L  J.  Morey,  V.  Vahedi,  3.  P.  Verboncoeur.  to  be  presented  at  this  con¬ 
ference. 

’M.  A  Lieberman,  3.  Appl  Phys.  85,  4186  (1989). 
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WinGraphics:  AN  OPTIMIZED  WINDOWING  ENVIRONMENT 
FOR  INTERACTIVE  REAL-TIME  SIMULATIONS 

John  P.  Vcrboncocur  and  Vahid  Vahedi 
Plasma  Theory  and  Simulation  Group 
Electronics  Research  Laboratory 
University  of  California 
Berkeley,  CA  94720 

We  have  developed  a  customized  windowing  environment,  WinGraphics,  which  provides  particle 
simulation  codes  with  an  interactive  user  interface.  The  environment  supports  real-time  animation 
of  the  simulation,  displaying  multiple  diagnostics  as  they  evolve  in  time.  In  addition,  keyboard 
and  printer  (PostScript  and  dot  matrix)  support  is  prov.ded.  The  simulation  codes  are  structured 
as  shown  in  figure  1. 


FIGURE  1 .  Schematic  representation  of  the  interaction  between 
WinGraphics  and  the  physics  kernel. 


The  physics  kernel  is  portable  to  any  machine  supporting  standard  C  The  INTT  module  scans  the 
input  file  containing  the  physical  parameters  of  the  problem  and  initializes  the  diagnostic  windows. 
INTT  also  sets  up  memory  for  array  storage.  The  environment  provides  hooks  for  the  physics 
kernel  to  run  continuously  (there  is  no  time  limit  •  die  code  can  run  indefinitely)  or  step  through 
individual  timesteps.  The  screen  is  refreshed  each  timestep,  and  all  user  requests  are  processed 
by  the  WinGraphics  MANAGER.  When  the  simulation  is  in  the  running  state,  the  MANAGER 
is  invoked  at  each  time  step  (if  not  running  the  MANAGER  is  constantly  invoked)  to  check  the 
keyboard  buffer  for  the  user’s  requests  (messages).  Pending  messages  are  translated  and  dis¬ 
patched  by  die  MANAGER  until  the  user  requests  QUIT. 

The  WinGraphics  environment  provides  menus  (as  shown  in  figure  2),  messages,  prompts,  and 
dialog  boxes  as  well  as  keyboard  support  for  navigating  the  above.  These  facilities  allow  interactive 
control  of  both  the  simulation  and  the  output  of  the  simulation. 


Figure  2.  The  WinGraphics  menus  for  P0W1  [1]. 

WinGraphics  controls  the  video  display  by  managing  ea ch  window  as  an  object  with  the  following 
attributes: 

•  Labels  for  each  axis. 

•  The  coordinates  of  the  outer  window  frame  and  the  client  area. 

•  Scales  for  the  x-  and  y-axes.  These  quantities  are  useful  for  normalized  codes  where 
values  must  be  converted  to  some  units  before  plotting. 

•  Maximum  and  minimum  x  and  y  coordinates  of  the  plot. 

•  State  of  the  object  The  environment  checks  the  state  to  determine  if  the  window  is 
on/off,  if  it  was  previously  on,  the  type  of  plot(s)  the  window  contains  (semilog,  linear, 
or  scatter),  if  die  window  is  overlapping  other  windows,  and  if  it  is  currently/previously 
being  autorescaled1.  The  state  is  used  to  minimize  redrawing,  objects  which  have  not 
changed  since  the  last  screen  update  are  not  redrawn. 


1.  Autorescale:  a  process  of  scanning  the  x  and  y  arrays  for  the  maximum  and  minimum  val¬ 
ues.  An  autorescaled  plot  fills  the  entire  client  area  of  its  window,  and  causes  the  axis  labels  to 
be  recalculated  each  time  they  change. 


•  Linked  list  of  curves  to  plot  The  contents  of  the  linked  list  include  the  number  of  points, 
pointers  to  the  arrays,  and  the  color  of  the  curve.  This  enables  windows  to  display 
multiple  curves  using  the  same  axes. 

By  manipulating  the  above  characteristics  interactively,  windows  can  be  opened,  closed,  moved 
to  a  new  location  on  the  screen,  made  smaller  or  larger  in  two  dimensions,  and  overlapped.  In 
addition,  WinGraphics  provides  an  interactive  crosshair  with  which  the  user  can  point  to  any 
location  on  the  screen  and  display  the  coordinates  of  that  point  in  the  units  of  the  plot  containing 
die  point 

WinGraphics  provides  management  of  the  segmented  Intel  80x86  memory  by  allocating  only  the 
required  space  for  arrays  and  structures.  Near  arrays  (addressed  using  the  DS  segment  register 
and  a  16-bit  offset)  are  used  for  performance -critical  variables,  while  far  arrays  (addressed  with 
a  16-bit  segment  and  a  16-bit  address)  are  used  for  large  arrays  which  cannot  fit  in  the  64kB  near 
heap  of  the  80x86. 

Since  the  code  can  run  for  an  indefinite  number  of  timesteps,  arrays  which  accumulate  some 
physical  quantity  in  time  must  be  managed.  The  histories  are  combed  when  the  array  limit  is 
reached,  retaining  only  every  nth  value.  This  leads  to  loss  of  temporal  resolution  as  the  simulation 
is  advanced  in  time.  To  retain  local  temporal  resolution,  local  arrays  are  used  which  contain  data 
for  contiguous  timesteps. 

The  codes  currently  running  in  the  WinGraphics  environment  include  ESI  (Electrostatic  1 
Dimensional  Periodic  Plasma  simulation)  [2],  PDW1  (Plasma  Device  Workshop  1  Dimensional 
Cartesian  bounded  simulation)  [1],  and  PDC1  (the  cylindrical  counterpart  of  PDW1)  [3]. 

Figure  3  shows  a  sample  PostScript  output  of  WinGraphics  for  ESI  displaying  the  variation  of 
potential,  4>(x)  (left),  and  the  time  history  of  electrostatic  field  energy  (right)  for  Landau  damping 
with  8192  particles. 


Figure  3.  PostScript  output  of  WinGraphics  for  ESI . 


As  a  performance  comparison  of  the  codes,  we  use  ESI  running  a  Landau  damping  simulation 
with  16000  electrons.  The  results  of  the  comparison  are  shown  below.  The  Cray  values  are  given 
using  the  NMFECC  Cray  computers  using  the  vectorized  FORTRAN  version  of  the  code  due  to 
A.  B.  Langdon  (LLNL). 


computer 

clock  speed  (MHz) 

psec  per  particle  per  push 

IBM  AT 

6 

1038 

IBM  PS/2  Model  80  16 

183 

Dell  310 

20  +  memory  cache 

118 

Cray  I 

vectorized 

1.8 

Cray  II 

vectorized 

1.4 

From  this  information,  we  can  extrapolate  for  new  hardware.  For  example,  an  80386/80387 
running  at  33  MHz  with  sufficiently  fast  memory  should  perform  at  7 1  psec  per  particle  per  push, 
an  80486  at  33  MHz  should  perform  at  24  psec  per  particle  per  push,  and  an  80486  running  at  50 
MHz  should  result  in  performance  of  16  nsec  per  particle  per  push. 

Clearly  we  are  rapidly  approaching  the  performance  needed  to  run  two  dimensional  codes  on 
PC’s.  Note  that  as  the  number  of  dimensions  increases,  the  relative  computational  cost  of  the 
field  solve  increases,  so  one  cannot  directly  extrapolate  from  the  ID  panicle  pushes.  Given 
sufficient  memory,  it  is  estimated  that  a  minimum  speed  of  40  psec  per  particle  per  push  (25,000 
particles  per  second)  is  required  to  run  a  2D  simulation  in  real  time.  Memory  requirements  vary 
widely  depending  on  the  problem,  but  a  good  minimum  estimate  is  800k  bytes  (500k  for  particle 
quantities  of  25,000  particles,  120k  for  quantities  on  a  100x100  grid,  etc.) 
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BOUNDED  PLASMA  DEVICE  SIMULATION  WITH  PDW1, 
INCLUDING:  EXTERNAL  RLC  CIRCUIT, 

DC  AND  RF  DRIVE,  AND  COLLISIONAL  PROCESSES. 


Ian  J.  Morey,  John  P.  Verboncoeur,  and  Vahid  Vahedi 
Plasma  Theory  and  Simulation  Group 
Electronics  Research  Laboratory 
University  of  California 
Berkeley,  CA  94720 

PDW1  (Plasma  Device  Workshop  1-dimensional),  running  in  the  WinGraphics  environment  [1], 
is  a  versatile  many-particle  bounded  plasma  simulation  code  used  to  study  diverse  plasma  and 
boundary  conditions.  The  boundaries  are  connected  through  an  external  RLC  circuit  with  time 
dependent  voltage/current  sources.  Having  a  versatile  input  deck  is  an  efficient  way  of  specifying 
plasma  parameters  to  the  code  without  recompiling.  In  PDW1,  input  deck  parameters  are  given 
in  MKS  units.  The  list  includes: 


plasma  variables: 

device  variables: 
circuit  variables: 
electrode  variables: 
collision  variables: 
numerical  variables: 


number  of  species,  characteristics  of  each  species  and  back¬ 
ground  charge  and  current  densities 
electrode  separation,  electrode  area  and  background  e 
R,  L,  C,  V(t)  and  /(/) 

absorption,  secondary  emission  and  plasma  sources 
neutral  pressure,  elastic,  excitation  and  ionization  collisions 
A r.  Ax,  number  of  cells,  charge  per  computer  particle  and  FFT 
period. 


PDW1  employs  four  distinct  circuit  solvers  to  handle  the  full  range  of  external  circuit  parameters. 
For  the  general  series  RLC  circuit  with  voltage  source,  a  2nd  order  backward  Euler  is  used.  For 
C  -»  0,  the  external  circuit  becomes  an  open  circuit,  and  external  parameters  are  ignored.  For 
C  and  R=L=0,  external  potentials  are  applied  directly  across  the  plasma  region.  The  final 
case  is  an  ideal  current  source,  which  drives  the  specified  current  from  one  plate  to  the  other 
(external  circuit  parameters  do  not  affect  an  ideal  current  source). 


Figure  1.  Time  and  frequency  responses  of  voltage  and 
current  signals  in  an  RF  discharge. 


FFT  routines  are  used  to  obtain  the  frequency  response  of  any  signal  accumulated  in  time.  The 
quantities  currently  transformed  include  the  driven  wall  potential,  the  potential  at  the  middle  of 
the  system  and  the  current.  The  total  impedance,  Z(f),  of  the  plasma  region  is  calculated  in  the 
frequency  domain  as  the  redo  of  V(f)  to  1(f).  An  example  is  shown  in  figure  1 .  for  an  RF  discharge. 

Collisions  with  neutral  particles  have  been  included  in  PDW1  so  partially  ionized  discharges  can 
be  The  present  model  can  be  used  for  elastic,  inelastic  and  ionizing  electron-neutral 

collisions,  and  can  be  extended  to  charge-exchange  ion-neutral  collisions.  The  full  3-D  character 
of  a  collision  is  modelled  with  just  two  velocity  components,  parallel  and  perpendicular  to  the 
spatial  dimension  of  the  simulation. 

The  neutrals  are  assumed  to  be  uniformly  distributed  between  the  boundaries  with  a  constant 
density.  The  probability,  Pm ,  that  the  m-th  electron  has  a  collision  with  a  neutral  is  given  by 

Pm  «  1  -expC-no^EJv^Ar)  (1) 


where  n  is  neutral  density,  Or  is  the  total  cross-section  (the  sum  of  cross-sections  for  elastic, 
excitation  and  ionizing  collisions)  and  v„  is  the  velocity  of  the  m-th  electron.  Laboratory 
cross-sections  are  entered  into  arrays  at  the  start  of  the  run.  A  collision  occurs  if  a  random  number 
R(  R  e  [ 0,1  D  is  less  than  P*.  Another  random  number  is  then  used  to  determine  what  type  of 
collision  has  occurred,  based  upon  the  following: 

R  £  a^Eya-riE)  '.elastic 

a+CEyOjfE)  <R  £ (o+(E)  +  am(E)yctt(E)  : excitation 

(o+(E) + G^(£))/Or(£)  <  R  tionization 


Once  the  type  of  collision  is  determined,  the  energy  of  the  scattered  electron  is  obtained.  Energy 
is  unchanged  for  elastic  collisions,  while  excitation  energy  is  lost  during  excitation  collisions. 
During  an  ionizing  collision  the  energy  is  divided  between  the  two  electrons  using  a  differential 
ionization  cross-section  of  the  form: 


S(E,T ) 


ME) 

T2+B\E) 


(2) 


where  £  is  the  incident  electron  energy,  7  is  the  scattered  electron  energy,  and  A  and  B  are  general 
functions  of  £.  This  is  a  simplified  form  of  the  cross-section  used  by  Peterson  and  Allen  [2]. 

Lastly,  the  directions  of  the  scattered  and  ejected  electrons  are  obtained  by  using  the  same  dif¬ 
ferential  cross-section  to  determine  the  scattering  angle  after  all  collisions,  which  has  been  taken 
from  Den  Hartog  et  al  [3],  and  has  the  form 

0(7,  x) _ 7 _  (3) 

0(7)  *  4tt(l  +  7  sin2(x/2)) ln(l  +  7) 

where  x  is  the  angle  between  the  incident  and  scattered  velocities. 


Figure  2.  RF  discharge  ion  and  electron  phase  space  (left)  and  potential  #*)  (right)  over 
one  half-cycle  with  right  wall  driven  at  1 0  MHz,  Ar * in  sec 


Figure  2  shows  a  sample  run  of  PDW1  for  an  RF-driven  discharge  over  one  half-cycle,  displaying 
die  velocity  phase  space  (velocity  versus  position  for  each  particle),  and  the  variation  of  potential, 
4(rX  from  wall  to  wall. 

The  left-hand  plot  in  figure  3  shows  the  velocity  phase  space  in  expanded  scale,  as  the  ions 
accelerate  out  to  the  walls  from  the  bulk  plasma  due  to  the  presence  of  an  average  positive  potential, 
as  seen  in  die  right-hand  plot,  in  the  middle  of  the  system. 


Figure  3.  Enlarged  to  show  ion  velocities  to  die  left;  time  response  of  mid-potential  to  the  right. 
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The  Traveling- Wave-Tube  Code  EBC 

by  LJ.  Morey  and  C.K.  Birdsall 
ERL,  University  of  California,  Berkeley,  CA  94720 


IBC  (Interactive  Beam-Circuit)  is  a  one-dimensional  many  particle,  traveling-wave-tube 
simulation  code  which  has  been  developed  to  run  interactively  on  a  PC  or  Workstation. 
The  code  follows  all  of  the  particles  in  the  tube,  rather  than  just  those  in  one  wavelength. 
This  allows  for  nonperiodic  inputs,  nonuniform  line  and  a  large  set  of  spatial  diagnostics 
which  can  be  displayed  graphically. 

The  code  uses  particle-in-cell  techniques  (Birdsall  and  Langdon,  1985)  to  model  the 
motion  of  beam  electrons  and  simple  finite  difference  methods  to  model  the  fields  of  the 
coupled  slow- wave  transmission  line.  The  coupling  between  the  beam  and  the  transmission 
line  is  based  upon  the  finite  difference  equations  of  Brillouin  (1949),  and  particle-in-cell 
techniques  are  also  used  to  include  the  space-charge  effects,  in  a  manner  similar  to  that 
used  by  Hess  (1961). 

The  transmission  line  equations  used  in  the  code  are: 


YL  -  vi~" 

A  t 

r*+l/2  rt— 1/2 

•*«-!/ a  ~~  *«- 1/2 
A  t 


1  Ij+i/2  ~~  Jj-1/2  ,  *(fc)  Qi,i  -  Qfci1 

C  A*  C  At 


iv?-vu 

L  Ax 


(1) 

(2) 


The  variables  V/  and  are  the  transmission  line  voltage  and  current  at  the  points  i  Ax 
and  (*+1/2)  Ax  respectively,  and  at  times  tAt  and  (<+l/2)A<  respectively.  The  capacitance 
and  inductance  per  unit  length  are  C  and  L.  The  beam  charge  within  a  section  of  the 
tube,  Qh,i,  induces  a  charge  on  the  transmission  line  «(&)<?*,,,  where  k  is  the  wavenumber 
of  the  space-charge  waves.  The  factor  «(i)  reflects  the  nature  of  the  space-charge  coupling 
constants  for  finite  diameter  TWT  and  depends  on  the  relative  diameters  of  the  beam  and 
the  waveguide,  and  the  transverse  field  profiles  of  the  different  wave  modes. 

From  the  equation  of  continuity,  another  driving  term  could  be  used,  similar  to  Pierce 

(1950): 

_  *(fc)  ~  1/2  ,  . 

C  Ax 

This  option  is  available  in  the  code.  Using  a  method  of  assigning  the  charge  of  each  electron 
particle  to  the  two  nearest  cell  boundaries  and  the  current  of  each  electron  particle  to  the 
two  nearest  cell  centers,  these  two  driving  terms  produce  the  same  growth  (satisfying 
continuity). 
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We  assume  that  the  space-charge  field  variables  have  a  radial  profile  of  the  form  J o (kxr) 
(TMoi  mode),  where  kx  =  2.405/i^m  and  R^m  is  the  beam  radius.  This  implies  that 
there  is  a  conductor  at  the  edge  of  the  beam  and  allows  Poisson’s  equation  for  the  space- 
charge  potential,  5,  and  density,  p,  in  cylindrical  coordinates  to  be  simplified  to 


P 

eo 


(4) 


where  ^(x,  r,  0)  =  4(z)Jo(kxr)  and  similarly  for  p.  We  also  assume  that  there  is  no  DC 
space  charge  because  the  unperturbed  electron  beam  is  neutralized  by  cold  immobile  ions. 
Therefore,  even  though  the  AC  space-charge  density  has  a  Jo(kxr)  profile,  we  can  choose 
a  uniform  radial  density  profile  for  the  beam  and  the  ions. 

An  additional  feature  has  also  been  added  so  that  the  effects  of  space  charge  can  be 
observed  easily.  The  constant  SC  coup  allows  the  user  to  change  the  strength  of  the  space- 
charge  force.  When  SCcoup  —  1,  the  full  effect  of  the  space  charge  is  included,  but  when 
SC  coup  =  0  the  space-charge  force  has  no  effect,  so 


E(total)i+ 1/2  —  E (ci'rcwtt) «-f 1  /2  +  SCcOUp  *  E[tpaee-charg€)i+l/2-  (5) 

The  same  relationship  applies  to  the  total,  circuit  and  space-charge  potentials,  and  the 
electron  bunching  observed  in  IBC  with  SCcoup  =  1  agrees  with  the  results  of  Hess 
(1961). 

The  transmission  line  is  terminated  at  both  ends  with  a  resistor  and  a  voltage  source. 
However,  these  terminations  can  be  matched  only  when  there  is  no  beam-circuit  coupling. 
The  greater  the  coupling  between  the  beam  and  the  circuit,  the  greater  the  reflected  signal. 
U  the  loop  gain  exceeds  unity,  the  tube  will  become  an  oscillator.  To  avoid  this,  a  section 
of  the  transmission  line  can  be  made  lossy,  or,  an  absorbing  termination  length  may  be 
added  to  the  tube  -  a  region  where  the  loss  increases  exponentially  and  the  beam-circuit 
coupling  decays  to  zero. 

The  boundary  conditions  for  the  particles  are  absorption  at  both  the  source  and  collec¬ 
tor.  The  beam  particles  are  injected  uniformly  at  the  source.  The  space-charge  potential 
is  fixed  at  zero  at  both  the  source  and  collector. 

When  the  code  is  run  on  an  IBM  PS2,  with  100-200  cells  and  up  to  5  particles  per 
ceH,  many  different  features  of  a  TWT  interaction  can  be  observed  in  a  one  hour  session. 
Quantities  such  as:  the  beam  density  and  current,  the  two  driving  terms,  the  circuit  voltage 
and  current  and  the  power,  can  be  displayed  individually  or  all  together.  The  following 
diagrams  show  the  output  for  a  tube  of  length  0.2  meters,  a  gain  of  C  =  0.056  and  an 
absorbing  termination. 
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PDC1:  ONE  DIMENSIONAL  RADIAL  CODE  FOR  A  CYLINDRICAL 
PLASMA  DEVICE  WITH  AN  EXTERNAL  RLC  CIRCUIT 


M.  Virginia  Alves1,  Vahid  Vahedi,  Charles  K.  Biidsall 
Plasma  Theory  and  Simulation  Group 
Electronics  Research  Laboratory 
University  of  California 
Berkeley.CA  94720 

PDC1  is  a  one  dimensional  (radial),  electrostatic  particle  simulation  code  with  the  plasma 
bounded  by  two  concentric  cylinders.  These  electrodes  can  be  coupled  to  an  external  circuit  as 
shown  in  Figure  1.  Only  spatial  variations  in  r  are  considered.  The  principles  applied  in  this 
code  are  the  same  as  those  used  in  the  planar  model  code  PDW1,  widely  used  since  1983  [1]. 


Figure  1.  Cylindrical  plasma  device  with  an  external  RLC  circuit. 

PDC1  is  very  flexible  and  runs  in  the  WinGraphics  environment  [2].  It  is  being  applied,  for 
example,  to  RF  discharges,  in  which  the  RF  powered  electrode  and  the  grounded  electrode  have 
different  areas,  and  to  Plasma  Immersion  Ion  Implantation.  For  convenience,  the  inner  cylinder 
is  considered  powered  and  the  outer  cylinder  is  grounded.  In  addition,  in  1D-3V  (r,  v„  v*  v,) 
with  an  applied  axial  magnetic  field,  B^,  many  magnetized  non-neutral  plasma  problems  may 
be  addressed,  with  or  without  the  center  cylinder. 

We  obtain  the  finite  difference  equations  for  fields  (£)  and  potential  (<}>)  related  to  density  (p), 
on  a  spatial  grid,  using  Gauss’  law,  in  order  to  guarantee  conservation  of  flux.  The  particles  are 
considered  to  be  cylindrical  shells,  uniform  along  z.  The  grid  quantities  are  indexed  by  j.  The 

•  '0 

grid  points  are  spaced  uniformly  in  r,  Ar  *  -jg~,  where  NC  is  the  number  of  cells,  and  r0  and 

r^c  are  the  inner  and  outer  cylinder  radius.  The  charge  weighting  guarantees  charge  conserva¬ 
tion. 
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Applied  to  the  cylindrical  surface  at  j+\.  Gauss’  law  produces  the  radial  electric  field  from  the 
charges  as  in  Birdsall  and  Langdon  [3]: 


—  =  2nr  XE  ,-2 lor  ,£  ,  where  r .  — -  and  j> 0 


(1) 


For  j  -  0  we  have: 


Q0 

—  =  2jtr,£,  -  lizrJEt 
e  si 


In  order  to  obtain  an  equation  in  we  use  single  cell  differencing:  * 

4>>-i— 4»y 


Ar 


'♦5 


where  Ar  ,  sr  ,-r 


Q.  , 

The  charge  density  is  obtained  from  py  =  where  Vy  s  tt(Ar)y  with: 

(Ar2)  =r2  ,-r2  ,  j>0 

1  J2  }~1 

(Ar2)0  =  r,2-r2 


(2) 


(3) 


(4) 

(5) 


Using  (3M5)  in  (1)  and  (2)  we  obtain  a  three-point  finite  difference  form  of  Poisson’s  equation: 

=  0.5)6  .-2^6+0: +0.5)6:.,  i> 0  (6) 


-^Ar2T  =  (0 -0.5»>>.l-2r^  +  (r>+0.5)(J)>+1  j  >  0 

€<► 


(r0+0.25)Ary+r0CT 
(r0+OJ)e 


(7) 


In  (6)  ry  were  normalized  by  A r,  and  to  derive  (7),  we  assume  £o  =  “,  where  a  is  the  charge 

surface  density  at  the  inner  cylinder.  Equation  (7)  is  one  of  the  boundary  conditions  used  to 
solve  the  Poisson  equation.  The  second  boundary  condition  is  obtained  assuming  that  die  outer 
electrode  is  at  the  reference  potential,  =0. 

Equations  (6)-(7)  provide  a  tridiagonal  matrix  which  is  solved  to  find  the  potentials  $y.  Then  the 
electric  fields  are  obtained  using: 

5—*'^ for  y  =  l,2,...,JVC-l 

The  electric  field  at  the  first  grid  point  is  obtained  from  (2);  for  the  last  grid  point  we  use  a 
similar  equation  for  EK. 


In  order  to  solve  equations  (6)-(7)  for  the  potentials,  the  charge  density  must  be  known  on  the 
same  grid.  The  weighting  used  to  accumulate  charge  is  similar  to  linear  weighting.  We  use  the 
area  of  tings  to  weight  the  charges  ([3];  without  variations  in  6). 

The  weighting  of  the  electric  field  to  the  particles  is  done  in  the  same  manner.  The  charge  accu¬ 
mulated  on  the  grid  point  at  the  boundary  is  some  fraction  (in  PDW1,  exactly  [1])  of  that 

which  would  be  accumulated  at  the  grid  point  which  is  not  at  a  boundary  if  the  physical  charge 
density  were  die  same  at  the  two  points.  To  compensate  for  this  we  use  a  volume  corresponding 
to  half  a  cell  (see  eq.(5))  when  we  define  charge  density  at  die  last  grid  points  (j=0,  j=NQ. 


In  order  to  load  the  system  at  t=0  with  a  uniform  density,  we  assume  a  constant  density  in  r 
equal  to  ^  w^cre  W  is  the  total  number  of  particles  to  load  the  system.  Since  the  density 

is  uniform  in  r,  at  each  position  r,  we  have: 


r i  -'’.-i  =' 


r 2  -2 
rNC  ~0 


where  r-.,  is  the  position  of  the  previous  particle.  The  algorithm  (8)  is  used  to  determine  the 
positions  of  the  particles  at  t~0. 


The  cylindrical  plasma  device  with  an  external  circuit  RLC  is  showed  in  Figure  1.  In  order  to 
solve  the  circuit  equation  we  use  die  2nd  order  backward  Euler  method  and  we  treat  four  differ¬ 
ent  cases  (as  in  PDW1  [4]): 


a)  C=0,  open  circuit; 

b)  current  source; 

c)  R  *  L  *  0  and  C  — »  short  circuit; 

d)  general  case:  RLC  circuit  with  voltage  source. 


It  is  important  to  consider  carefully  the  relation  between  o  and  ^  in  order  to  specify  the  correct 
boundary  condition  for  solving  Poisson’s  equation. 

Cylindrical  computer  experiments  were  performed  many  years  ago  with  ID  models  in  order  to 
ascertain  the  possibility  of  electrostatic  inertial  confinement  [5]  and  the  dynamics  of  limiting 
charge  and  current  [6]. 

The  following  plots  illustrate  a  run  of  PDC1  for  an  RF  discharge  driven  by  a  current  source 
displaying  the  velocity  phase  space  (velocity  versus  position  for  each  particle),  die  potential 
across  the  bounded  region  as  a  function  of  position,  and  the  time  history  of  the  potential  at  the 
inner  electrode. 
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Multi-Scale  Particle  Simulation  of  Bounded 

Plasmas 
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1  Introduction 

We  are  using  the  multi-scale  technique  [l]  to  model  bounded  systems.  Certain  bounded  systems 
are  a  naturally  suited  for  the  multi-scale  method  because  the  boundary  layer  (or  sheath)  that 
forms  at  the  wall,  which  is  usually  a  short  spatial  (~  A  #e)  and  time  scale  (~  Wp*)  structure, 
can  significantly  affect  the  bulk  plasma  behavior.  One  goal  is  to  understand  the  interaction 
between  the  bulk  plasma  and  the  sheath.  If  the  relevant  short  time  scale  physics  is  local  to  a 
few  (or  one)  known  spatial  regions,  then  one  can  take  advantage  of  this  by  advancing  particles 
with  variable  At  depending  on  position,  hence  reducing  computing  time.  The  unmagnetized 
sheath  problem  is  such  a  case. 

The  model  is  a  one  dimensional  bounded  slab  with  kinetic  ions  and  electrons.  We  start  with 
a  collisionless  and  unmagnetized  system  for  simplicity.  The  right  boundary  is  a  conducting  wall 
that  absorbs  all  particles  that  come  in  contact  with  it.  The  left  boundary  is  a  symmetry  point, 
where  the  particles  are  reflected.  We  allow  a  specified  initial  distribution:  f(x,  v,t  =  0). 

In  order  to  test  the  numerics  of  both  the  multi-scale  method  and  boundary  conditions  we  are 
using  the  following  test  problem:  a  cutoff  Maxwellian  distribution  for  the  electrons  and  fixed  (or 
infinitely  massive)  ions.  The  system  has  an  analytic  solution,  so  the  run  may  be  started  from 
equilibrium.  This  gives  us  a  benchmark  and  tests  the  fast  time  scale  electron  sheath  dynamics. 
Results  using  variable  At  (with  ratio  of  smallest  to  largest  of  1:128)  will  be  given.  In  the  future, 
we  intend  to  use  the  more  general  model  to  study  time  dependent  bounded  plasma  problems, 
such  as  a  plasma  expanding  toward  a  conducting  wall. 

2  The  Multi-Scale  Method 

The  goal  is  to  model  bulk,  macroscopic  (u  <  up,  A  >  Xq)  plasma  phenomena  globally,  while 
also  accurately  treating  short  time  (~  w"1)  and  short  space  scales  (~  A*))  in  local  regions 
where  microscopic  physics  is  important  (e.g.  sheaths).  We  accomplish  this  by  allowing  groups 
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of  particles  to  move  at  different  At’s.  The  Direct  Implicit  method  is  used  to  allow  for  large  upAt 
[2].  Each  group  of  particles  (call  them  Gm),  is  pushed  every  2m  time  steps,  using  Atm  -  2 mSt, 

m  =  0,1, 2, _ ,  where  St  is  the  smallest  time  increment.  In  order  to  avoid  having  “special’’ 

time  steps,  we  do  not  move  all  the  particles  of  a  given  group  at  once.  We  define  subgroups 
of  Gm  called  blocks  (call  them  Blm).  Given  a  group  Gm,  there  are  2m  blocks  Blm,  where 
l  =  0,1,2, ...,  (2™  -  1).  Each  block  in  a  group  is  moved  at  a  different  time  step.  'OVe  move  a 
block  Blm  at  time  step  n,  with  t  =  n St,  when:  (n)mod(2m)  s  1.  As  an  example,  consider  3 
groups  (m  =  0,1,2): 


Group 

Step  size 

Block 

When  pushed 

Go 

Ate  =  St 

OO 

every  step 

G\ 

II 

** 

< 

B f 

even  steps 

B\ 

odd  steps 

G-i 

II 

<1 

B°2 

if  n  mod  4=0 

B\ 

if  n  mod  4=1 

Bi 

if  n  mod  4=2 

B\ 

if  n  mod  4=3 

To  illustrate  how  particles  axe  moved,  take  time  step  n=7: 


Three  blocks  are  moved  at  n=7:  B§,  B\  and  B\  (represented  by  the  solid  arrows).  Contributions 
to  the  free  streamed  charge  density  p,  and  the  “effective”  susceptibility  x  from  aD  blocks  are 
known  in  advance  of  the  step;  those  from  blocks  B$,B%,  B\  and  B\,  have  been  obtained  by 
interpolation  in  time.  Thus,  the  field  is  known  before  the  blocks:  Bjj,  B}  and  2?f  are  advanced 
to  n=7. 


2 


3  Physics  of  the  Test  Problem 


The  sheath  problem  we  use  is  a  simple  example  to  test  the  electron  sheath  dynamics.  f(x,  v) 
for  the  electrons  is  a  cut-off  Maxwellian: 


/(*,») 


^  tiqCq  exp  \V\<Ve(x), 

[  0;  M  >  ®e(x). 


(1) 


where  n0  =  ne(z  =  0),  mv\  —  T,  cq  —  ^2y/2vreri and,  vq  =  vc(x  =  0).  The  ions 

y/2  V  g(f> 

are  fixed:  n,(x)  =  no  Changing  variables;  u  =  — - ,  and  ifr  =  — ,  Poisson’s  equation  becomes: 

2  vj  T 


1  f  eriy/uo  +  t/> 
d^~  A^\  ertytiS 


We  can  now  solve  Eq.  (2)  with  the  appropriate  boundary  conditions,  and  knowing  <$>{x)  we 
load  the  particles  (electrons)  according  to  Eq.  (1). 


4  Comparison  of  Computing  Time 

Comparison  of  cpu  time  for  two  runs:  one  with  1  At  group  and  the  other  with  8  At  groups.  No. 
of  particles  =  100,000;  no.  of  timesteps  =  10,000.  Time  is  given  in  microseconds  per  particle 
per  times  tep. 


largest  At  . 
smallest  At  ' 

1 

128 

Gain 

total  time : 

4.62 

1.56 

2.96 

without  io : 

3.75 

0.670 

5.60 

without  io,  or  initialization.  : 

3.63 

0.546 

6.64 

We  have  not  yet  fully  vectorized  the  code.  There  is  still  room  for  optimization,  although  the 
logic  involved  in  changing  particles  from  block  to  block  makes  vectorization  more  complicated. 
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Symmetry  Plane 


Conducting  Wall 


Figure  1:  Schematic  of  phase  space  for  the  bounded  multiscale  simulation.  The  case  shown  has 
6  groups.  Particles  change  At  groups  as  they  move  in  phase  space.  The  solid  line  represents  the 
aeparatrix,  ec(x),  where  inside  the  particles  are  trapped  and  outside  the  particles  hit  the  wall 
and  are  absorbed.  The  ion  are  fixed  and  the  electrons  have  a  cut-off  Maxwellian  distribution. 
The  conducting  wall  is  left  floating,  V/. 
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1  Introduction 

We  propose  a  computational  method  to  self-consistently  model  small  angle  collisional  ef¬ 
fects.  This  method  may  be  added  to  standard  Particle-In-Cell  (PIC)  plasma  simulations 
to  include  collisions,  or  as  an  alternative  to  solving  the  Fokker-Planck  (FP)  equation  using 
finite  difference  methods.  The  distribution  function  is  represented  by  a  large  number  of 
particles.  The  particle  velocities  change  due  to  the  drag  force,  and  the  diffusion  in  velocity 
is  represented  by  a  random  process.  This  is  similar  to  previous  Monte-Carlo  methods  [1,2], 
except  we  calculate  the  drag  force  and  diffusion  tensor  self-consistently.  The  particles  are 
weighted  to  a  grid  in  velocity  space  and  the  associated  “Poisson  equations”  are  solved  for 
the  Rosenbluth  potentials.  The  motivation  is  to  avoid  the  very  time  consuming  method  of 
Coulomb  scattering  pair  by  pair.  First  the  approximation  for  small  angle  Coulomb  collisions 
is  discussed.  Next,  the  FP-PIC  collision  method  is  outlined.  Then  we  show  a  test  of  the 
particle  advance  modeling  an  electron  beam  scattering  off  a  fixed  ion  background. 


2  Small  angle  Coulomb  collisions 


The  FP  equation  for  describing  small  angle  Coulomb  collisions  can  be  solved  numerically 
using  finite  difference  techniques.  An  alternate  method  is  to  follow  the  evolution  of  a  large 
number  of  particles.  Using  the  notation  of  Trubnikov  [3],  an  infinitestimal  “cloud”  of  test 
particles  can  be  represented  by  the  quantities  <  Av;  >  and  <  A >  to  the  same  order 
of  accuracy  as  the  FP  equation[4j,  and  are  defined  by: 


<  A  Vi  >= 


dvj(t) 

dt 


<  Ad.Auj  >=  ^{(c,  -  vi){vj  -  Uj)} 


(1) 

(2) 


The  bars  represent  averages.  These  may  be  interpreted  as  ensemble  averages  over  many 
initial  states,  or  as  an  ordinary  average  by  letting  the  number  of  particles  in  a  cloud,  N. 
become  large.  We  follow  the  second  interpretation  allowing  a  simpler  simulation  method. 
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Si  is  now  the  average  velocity  of  the  local  cloud.  The  i  and  j  subscripts  are  the  3  velocity 
components  in  Cartesian  coordinates,  (t>j  =  vx,  v-i  —  vy,  =  vx).  <  Av,  >  is  the  drag 
force  felt  by  an  average  particle  in  the  cloud.  <  An, A Vj  >  is  the  spreading  or  diffusion  in 
velocity  space  of  the  cloud. 

Numerically  we  advance  the  velocity  of  the  particles  by  using  the  following  equation  [4]: 


A Vi{tn)  =<  A V{  >  At  +  BijAWi(tn) 


(3) 


where  Btj  satisfies. 


BtkBitj  —  ^  Ar.A v j  > 


(4) 


where  AHr,-(tn)  =  VT,(tn)  -  W,(tn_1),  and  W,(t)  is  a  vector  function  that  represents  a 
Wiener  process  (or  Brownian  motion)  having  the  following  properties  <  Wi(t)  >=  0  and 
<  Wj(<)Wj(f)  >=  t6ij.  B  and  W  are  not  unique  and  only  need  to  satisfy  the  above  criteria. 
In  our  model  we  choose: 

AW,(t")  =  VAtR,  (5) 


where  are  independent  random  numbers  satisfying  <  R,  >=  0,  and  <  R,Rj  >=  We 
have  found  a  a  matrix  B  such  that  BBr  =  A ,  given  by: 


By  letting  AXJ  =<  Av,Av}  >,  we  note  that  Eq.  (6)  satisfies  Eq.  (4)  and  can  therefore  use 
this  B  in  Eq.  (3). 

The  “field”  quantities  are  obtained  from  two  functions  <t>  and  the  Rosenbluth  poten¬ 
tials,  which  in  turn  are  solved  by  the  two  Poisson  equations  [3]: 


vy  =  />,) 


(T) 


vV  =  **(*.•) 


(8) 


where  0  is  the  superscript  representing  the  field  species,  0  =  (i,  e).  Then  <  A r<  >  and 
<  Av,Avj  >  are  obtained  in  terms  of  the  Rosenbluth  potentials  using  the  following  equation 
(see  reference  [3]  for  a  derivation): 


<  Av;  >°=  - 


TQ/3d<t>a 

dVi 


(9) 
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(10) 


<  Avi&Vj  >°=  -2^  La/0 
0 


d2'i>0 

dvidvj 


Lo/0  _  \  ^^0^ 

is  a  constant  (given  here  in  cgs  units),  and  A  is  the  Coulomb  logarithm. 
a  and  0  represent  the  test  and  field  particles,  respectively  (e.g.  a  =  (t,e)  and  0  =  (i,e)). 


3  Outline  of  the  Method 

For  the  particle  advance,  we  start  with  the  simplest  scheme  [1],  similar  to  Euler’s  method, 
except  there  is  an  added  diffusion  term  (the  last  term  on  the  right): 

r"+1  =  «?+  <  At’,  >n  At  +  B”}y/AtRj  (11) 

Rj  are  independent  random  numbers  having  the  two  properties  given  above.  Following  a 
methodology  similar  to  PIC  simulation,  we  weight  the  particles  to  a  grid,  calculate  the  field 
quantities,  advance  the  particles  and  repeat  the  process.  The  basic  algorithm  is  as  follows: 

Step  1.  Weight  the  particles  to  a  grid  in  velocity  space  using  linear  interpolation. 

Step  2.  Solve  Eqs.  (7)  and  (8)  on  the  grid  for  <f>  and  rj.\ 

Step  3.  Obtain  <  A r,  >n  and  on  the  grid  using  Eqs.  (6),  (9)  and  (10). 

Step  4.  For  each  particle,  obtain  <  Ac,-  >"  and  5”  by  interpolating  from  the  grid  to 
the  particle  location  u;.  Then,  update  the  velocity  using  Eq.  (11). 


4  Test  of  the  particle  advance 


As  an  initial  test  of  the  particle  advance,  Eq.  (11),  we  model  an  electron  beam  scattered 
by  infinitely  massive  cold  ions.  For  this  test  case  <  Au,  >  and  <  Av.Arj  >  are  given  by 
the  two  simple  analytic  expressions: 


<  Av,  >=  -C-l 
v3 

<A^>-c(&-^) 


(12) 

(13) 


_  n0L° 
4  ?r 


where  C  =  — ; - .  We  do  not  calculate  the  field  quantities  on  the  grid  (step  1  and  2 

above),  but  rather,  use  Eqs.  (12)  and  (13)  directly.  This  provides  a  good  test  of  the  particle 
advance,  Eq.  (11).  We  expect  the  total  x-momentum  to  be  given  by  [3]: 


Px,toial(t )  =  Px,toial(t  —  0)CX p  (  3  t\ 

\  v0  / 


(14) 
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where  vo  is  the  initial  beam  velocity.  We  have  made  a  test  run  with  the  following  parameters: 
C  =  1,  vo  =  1,  VTe  =  0.01  and  At  =  0.001.  The  energy  error  is  less  than  0.2  percent  and 
the  momentum  error  is  less  then  1  percent  for  a  total  time:  T=l. 

5  Future  Work 

We  are  developing  a  code  to  test  the  FP-PIC  method  outlined  above.  In  the  near  future  we 
plan  to  implement  the  algorithm  in  a  Id  3v  electrostatic  code  in  order  to  study  combined 
collective  electrostatic  and  collisional  effects.  With  this  code  we  would  assume  spatial 
homogeneity  for  calculating  the  collisional  terms.  If  this  assumption  is  not  valid  one  would 
have  to  partition  f(v)  into  spatial  regions  j,  and  calculate  / }(v)  for  each  spatial  region.  One 
would  have  to  ensure  that  enough  particles  in  each  spatial  zone  to  adequately  fill  out  the 
distribution  function. 

In  addition,  we  need  to  address  the  accuracy  and  stability  issues  associated  with  this 
method,  and  study  the  feasibility  of  using  a  more  accurate  particle  advancing  scheme  than 
that  given  by  Eq.  (11). 
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1  Introduction 

The  transport  and  confinement  of  charged  particles  in  an  auxiliary-heated  mirror  plasma  is 
modeled  via  the  bounce-averaged  Fokker-Planck  (F-P)  code  SMOKE  [2]  and  the  direct-implicit 
partide-in-cell  (PIC)  code  TESS  [1].  The  test  case  studied  is  that  of  a  tandem  mirror  end  plug 
which  is  heated  by  the  injection  of  ECRH  at  both  the  fundamental  and  second  harmonic  of  the 
electron  gyrofrequency  =  «<*  *nd  wrf,i  —  2u/c*).  Figure  1  shows  the  magnetic  field  and 

potential  profiles  that  are  prescribed  for  the  test  case.  Both  electron-electron  and  electron-ion 
collisions  are  included.  Each  code  employs  a  relativistic  description  of  the  electron  dynamics  in 
one  spatial  dimension.  Each  code  follows  the  system  to  equilibrium,  where  the  loss  of  particles 
(resulting  from  collisions!  detrapping  of  the  loss-cone  distribution)  balances  the  trapping  of 
particles  (resulting  from  the  preferential  increase  in  the  magnetic  moment  of  the  particle  due 
to  the  resonant  absorption  of  RF  wave  energy). 

The  modeling  comparison  is  divided  into  three  sections.  The  first  entails  benchmarking  the 
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Figure  1:  Magnetic  field  and  potential  profiles  in  the  end  plug  of  a  tandem  mirror  plasma.  The 
central  cell  is  to  the  left,  and  the  end  wall  is  to  the  right. 

physics  results  from  the  PIC  code  against  those  from  the  F-P  code.  The  test  case  determines 
the  confinement  and  transport  of  electrons  only.  The  magnetic  and  electrostatic  fields  are 
prescribed  quantities.  Comparison  is  made  of  the  electron  velocity-space  density  contours  at 
different  axial  locations,  as  well  as  the  electron  density  and  kinetic  energy  profiles.  In  order  to 
determine  the  effect  of  particle  counting  statistics  on  these  quantities,  the  PIC  code  is  run  with 
varying  numbers  of  simulation  electrons.  The  PIC  code  is  also  run  including  the  ions,  and  a 
self-consistent  calculation  of  the  electrostatic  potential.  This  last  run  allows  us  to  determine 
the  accuracy  of  the  assumed  potential  profile  used  in  the  other  runs.  Next,  a  computational 
cost  analysis  is  performed  on  each  code.  In  addition  to  comparing  the  total  run  time  for  each 
of  the  codes,  various  code  dependent  costs  axe  provided,  such  as  the  cost  to  push  a  particle 
each  time  step  in  the  PIC  code  and  the  cost  to  "invert”  the  phase-space  matrices  in  the  F-P 
code.  Finally,  the  advantages  and  disadvantages  of  each  code  are  discussed  in  the  context  of 
the  chosen  test  case,  including  a  comparison  of  the  effort  involved  in  setting  up  the  input  deck 
for  each  code. 

2  The  Multiregion  Bounce-Averaged  Fokker-Planck  Code 

The  SMOKE  code  determines  the  evolution  of  a  particle  distribution  function  /  resulting  from 
Coulomb  collisional  and  RF-induced  velocity  space  diffusion.  The  code  numerically  solves  the 
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relativistic  Fokker-Planck  kinetic  equation 

f +v  ^+^KE+7vxB)/+r“',+r^]=o'  (i) 

where  t  is  the  time,  x  is  the  position,  v  is  the  velocity,  p  =  7 mv  is  the  momentum,  7  = 
(1  -  v’/c’r'/j  is  the  relativistic  factor,  m  is  the  rest  mass,  c  is  the  speed  of  light,  E  and  B 
are  of  the  electric  and  magnetic  fields  (static  and  RF  components),  and  Ta,u  and  Trf  are  the 
fluxes  arising  from  Coulomb  collisional  and  RF  wave  interactions  respectively.  Since  the  bounce 
time  of  particles  in  the  mirror  (j*)  is  much  smaller  than  the  collisional  or  RF  diffusion  times 
(T0olltTrf),  a  bounce-averaged  version  of  the  F-P  equation  [3]  is  solved  in  (c,/i)  phase  space 

= 1{D“% + + DJ) +iiD-rc+ + D“<)  ■  m 

where  {  =  (7  —  ljmc2  +  q$  is  the  total  energy,  ft  =  p^/2mfio  is  the  magnetic  moment  and 
n  =  /  dzj |vi|  |  is  the  particle  bounce  time.  Each  of  the  coefficient  Da0  have  contributions  from 
Coulomb  collisional  and  RF  wave  interactions.  The  quasilinear  model  of  Bernstein  and  Baxter 
[3]  is  used  to  model  the  RF-induced  diffusion. 


The  bounce-averaged  F-P  equation  (2)  is  solved  for  the  various  regions  (trapped-particle 
populations)  via  a  two-step  process.  First,  the  (c,  n)  phase  space  is  mapped  into  a  rectangular 
(Cartesian  coordinate)  region  (z,  y).  Second,  a  Galerkin  finite-element  method  is  used  to  solve 
(2)  in  the  (z,  y)  phase  space.  (The  same  procedure  is  used  to  solve  Poisson’s  equation  for  the 
Rosenbluth  potentials.) 


3  The  Direct-Implicit  Particle-In-Cell  Code 


The  TESS  code  computes  the  trajectories  of  individual  charged  particles  in  either  a  prescribed 

or  self-consistent  electrostatic  potential,  and  a  prescribed  magnetic  field.  The  implicit  nature  of 

the  code  allows  one  to  simulate  long-wavelength,  low-frequency  phenomena,  without  having  to 

resolve  high-frequency  effects,  such  as  electron  plasma  oscillations.  A  relativistic  guiding  center 

formulation  of  the  equations  of  motion  in  one  spatial  dimension  (£)  is  used 

dz  p, 

dt  -jm 


dpt 

dt 


qE,  -  ?-VBt  + 

7 


dp*\  dps\ 

dt  /  eoll  dt  '  rj 


dp  _  dp\  ,  dp\ 
dt  dt'eall  dt/rf 


(3) 
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In  general,  the  code  simulates  the  transport  and  confinement  of  both  ions  and  electrons.  The 
self-consistent  ambipolar  potential  is  then  calculated  via  the  direct- implicit  form  of  Poisson’s 
equation 

—  V  ‘  (l  +  !>•)  =  4*e  ( Zn «(*)  -  **«(*)) ,  (4) 

where  the  summation  is  over  species  *  (ions  and  electrons),  ns  is  the  free-streaming  or  explicit 
density  of  species  s  and  is  the  implicit  susceptibility  which  accounts  for  the  implicit  correction 
to  the  free-streaming  density.  The  equations  of  motion  (3)  and  modified  Poisson  equation  (4) 
are  solved  via  standard  finite  difference  techniques.  A  self-consistent,  relativistic,  Monte  Carlo 
binary  collision  model  [4]  and  a  quasilinear  RF  diffusion  model  [5]  (after  Bernstein  and  Baxter 
[3],  as  implemented  by  Rognlien  [6])  are  also  included. 
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